#Alexander F. Gazmararian
#afg2@princeton.edu
#January 9, 2024

#Purpose: Estimate primary models, robustness tests, and moerating effect of visibility

#Note: The code in this file needs to be sequentially run.

# Load packages
library(MatchIt)
library(tidyverse)
library(fect)
library(fixest)
library(tidylog)
library(sensemakr)
library(CBPS)
library(modelsummary)
library(panelView)
library(cobalt)
library(interflex)
library(kableExtra)
library(lmtest)
library(sandwich)
library(marginaleffects)
library(broom)
library(DirectEffects)
library(plm)
library(scales)
library(viridis)
library(gridExtra)
library(lfe)
library(here)
library(patchwork)

#Function to format latex tables
source(here("code", "fun", "fix_txt.R"))
#Function to plot the dynamic treatment effects alongside the equivalence test
source(here("code", "fun", "plot_fect.R"))
#Function for handling rounding when creating balance tables
source(here("code", "fun", "digitsByRows.R"))
#Function to create table output from DiD estimates via the FEct method
source(here("code", "fun", "create_fect_tab.R"))

# Set seed for replication
set.seed(10)

# Load data
g <- readRDS("data/output/pres_out.rds")
#Load matched data
g.nn <- readRDS("data/output/matched_df.rds")

# Define model for covariate balancing
f.match <- y ~ white + hispanic + foreign + college + income99pclog + poverty + rural_nonfarm + poplog + under40 + work_female

# Table set up
gof_map[gof_map$raw=="nobs",]$clean <- "$N$"
gof_map[gof_map$raw=="adj.r.squared",]$clean <- "Adjusted $R^2$"

# Define covariate names for creating tables
coefmap <- c(
  "coal_bin"="Coal",
  "post_shock"="Post-Shale Shock",
  "treated"="Coal $\\times$ Post-Shale Shock",
  "oilgas_prop"="Hydraulic Fracturing Employment",
  "unionshare"="Coal Union Share",
  "post_shock:unionshare"="Coal Union Share $\\times$ Coal Union Share",
  "white_s"="White",
  "post_shock:white_s"="White $\\times$ Post-Shale Shock",
  "hispanic_s"="Hispanic",
  "post_shock:hispanic_s"="Hispanic $\\times$ Post-Shale Shock",
  "foreign_s"="Foreign-born",
  "post_shock:foreign_s"="Foreign-born $\\times$ Post-Shale Shock",
  "college_s"="College",
  "post_shock:college_s"="College $\\times$ Post-Shale Shock",
  "income99pclog_s"="Income per capita (log)",
  "post_shock:income99pclog_s"="Income per capita (log) $\\times$ Post-Shale Shock",
  "poverty_s"="Poverty",
  "post_shock:poverty_s"="Poverty $\\times$ Post-Shale Shock",
  "rural_nonfarm_s"="Rural",
  "post_shock:rural_nonfarm_s"="Rural $\\times$ Post-Shale Shock",
  "poplog_s"="Population (log)",
  "post_shock:poplog_s"="Population (log) $\\times$ Post-Shale Shock",
  "under40_s"="Under 40 years",
  "post_shock:under40_s"="Under 40 years $\\times$ Post-Shale Shock",
  "work_female_s"="Female workforce",
  "post_shock:work_female_s"="Female workforce $\\times$ Post-Shale Shock",
  "debt06"="Debt-to-Income",
  "post_shock:debt06"="Debt-to-Income $\\times$ Post-Shale Shock"
)

# Coal country electoral reversal ----
# Claim on p. 1
# "In 2000, coal mining counties were 2.5\% more Democratic on average than the rest of the nation,
#casting their ballots for Al Gore, a Democrat who promised to tackle climate change.
#Since then, coal mines have shuttered. In 2020, coal mining counties were 10.3\% more Republican on average
#than the rest of the country, voting for Donald Trump, a Republican who pledged to roll back environmental regulations."
g %>%
  filter(year==2000 | year == 2020) %>%
  group_by(adjgroup, year) %>%
  summarize(gop = mean(RepVotesMajorPercent, na.rm = TRUE)) %>%
  filter(is.na(adjgroup) | adjgroup == "Coal") %>%
  mutate(adjgroup = ifelse(is.na(adjgroup), "nation", "coal")) %>%
  pivot_wider(names_from = c(adjgroup, year), values_from = gop) %>%
  mutate(
    gap_2000 = coal_2000 - nation_2000,
    gap_2020 = coal_2020 - nation_2020
  )

# Fig. 2: Pre-trends Prior to Matching----
# Calculate sample size
g %>% filter(!is.na(coal_dd) & !is.na(RepVotesMajorPercent)) %>% nrow()
# Create Figure 2 plot
p.line <- g %>%
  group_by(year, coal_dd) %>%
  summarise(gop = mean(RepVotesMajorPercent,na.rm=TRUE)) %>%
  mutate(coal = case_when(coal_dd == 1 ~ "Coal", coal_dd == 0 ~ "No Coal")) %>%
  ggplot(aes(x = year, y = gop, color = coal)) +
  geom_line(aes(lty = coal), linewidth = 1.25) +
  geom_point(aes(shape = coal), size = 3) +
  scale_color_grey() +
  scale_x_continuous(breaks = seq(1968, 2020, 4)) +
  scale_y_continuous(limits = c(30, 80), labels = percent_format(scale=1)) +
  annotate("text", x = 2018, y = 78, label = "Coal", size = 6) +
  annotate("text", x = 2018, y = 63, label = "No Coal", size = 6, color = "grey40") +
  labs(x = "Election", y = "Two-Party Republican Vote Share", shape = "", lty = "", color = "") +
  theme_bw(base_size = 14) +
  theme(
    panel.grid = element_blank(),
    legend.position = ""
  )
ggsave(
  p.line,
  device = "tiff",
  file = here("output", "figures", "fig2_pretrends_prematch.tiff"),
  width = 6.5,
  height = 3, 
  dpi = 300,
  scale = 1.5
)
ggsave(
  p.line,
  file = here("output", "figures", "fig2_pretrends_prematch.png"),
  width = 6.5,
  height = 3, 
  dpi = 300,
  scale = 1.5
)

# Main DiD Estimates ----
# Estimate using FEct
m.basic <- fect(RepVotesMajorPercent ~ treated + oilgas_prop,
                data = g.nn,
                index = c("county_state", "year"),
                force = "two-way",
                method = "fe",
                CV = TRUE,
                se = TRUE,
                nboots = 2000,
                parallel = TRUE,
                na.rm = FALSE)
m.basic
# Estimate placebo test
m.basic.placebo <- fect(RepVotesMajorPercent ~ treated + oilgas_prop,
                        data = g.nn,
                        index = c("county_state", "year"),
                        force = "two-way",
                        method = "fe",
                        CV = TRUE,
                        se = TRUE,
                        placeboTest = TRUE,
                        placebo.period = c(-2,0),
                        nboots = 2000,
                        parallel = TRUE,
                        na.rm = FALSE)
m.basic.placebo
## Fig. 3 ----
# Create dynamic treatment effects plot next to equivalence test
p.out <- plot_fect(fect.model = m.basic, equiv.p = "")
ggsave(
  p.out,
  device = "tiff",
  filename = here("output", "figures", "fig3_fect_main.tiff"),
  width = 6.5,
  height = 2.8,
  dpi = 300,
  scale = 1.5
)
ggsave(
  p.out,
  filename = here("output", "figures", "fig3_fect_main.png"),
  width = 6.5,
  height = 2.8,
  dpi = 300,
  scale = 1.5
)
## Table D3 ----
create_fect_tab(m.basic, m.basic.placebo, "/si_tab_D3_fect.txt")

# DiD Robustness ----
## Matrix Completion Estimator ----
#Estimate DiD estimates using MC estimator
m.mc <- fect(RepVotesMajorPercent ~ treated + oilgas_prop,
             data = g.nn,
             index = c("county_state", "year"),
             force = "two-way",
             method = "both",
             CV = TRUE,
             se = TRUE,
             nboots = 2000,
             parallel = TRUE,
             na.rm = FALSE)
# Estimate placebo using MC estimator
m.mc.placebo <- fect(RepVotesMajorPercent ~ treated + oilgas_prop,
                     data = g.nn,
                     index = c("county_state", "year"),
                     force = "two-way",
                     method = "mc",
                     placeboTest = TRUE,
                     lambda = m.mc$lambda.cv, 
                     CV = 0, 
                     se = TRUE,
                     placebo.period = c(-2,0),
                     nboots = 2000,
                     parallel = TRUE,
                     na.rm = TRUE)
### Table D.4 ----
create_fect_tab(m.mc, m.mc.placebo, "/si_tab_D4_fect_mc.txt")
### Fig. D.3 ----
m.mc$test.out$tost.equiv.p # Check p-value
p.mc.out <- plot_fect(fect.model = m.mc, equiv.p = "p<0.001")
ggsave(
  p.mc.out,
  filename = here("output", "figures", "si_fig_D3_fect_mc.png"),
  width = 6.5,
  height = 2.8,
  dpi = 300,
  scale = 1.5
)

## Sensitivity Analysis of DiD Estimates ----
### Table D.5 ----
#Linear Regression of Two-Party Republican Presidential Vote Share on the Interaction of Coal County Treatment Group Membership and Post-Shale Shock, 1972--2020
# Estimate cross-sectional model since covariates needed for benchmarking
m.ols <-
  lm(
    formula = RepVotesMajorPercent ~ treated + coal_bin + post_shock + factor(year) + state +
      + oilgas_prop +
      post_shock * (white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s +
                      rural_nonfarm_s + poplog_s + under40_s + work_female_s),
    data = g.nn
  )
m.ols.cl <- coeftest(m.ols, vcov. = vcovCL(m.ols, type = "HC1", cluster = ~ county_state))
# Add control for unionization
m.ols.union <- lm(
  formula = RepVotesMajorPercent ~ treated + coal_bin + post_shock + factor(year) + state +
    + oilgas_prop + unionshare * post_shock +
    post_shock * (white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s +
                    rural_nonfarm_s + poplog_s + under40_s + work_female_s),
  data = g.nn
)
# Add control for global financial crisis exposure
m.ols.gfc <- lm(
  formula = RepVotesMajorPercent ~ treated + coal_bin + post_shock + factor(year) + state +
    + oilgas_prop + debt06 * post_shock +
    post_shock * (white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s +
                    rural_nonfarm_s + poplog_s + under40_s + work_female_s),
  data = g.nn
)
# Estimate model with all controls
m.ols.all <- lm(
  formula = RepVotesMajorPercent ~ treated + coal_bin + post_shock + factor(year) + state +
    + oilgas_prop + unionshare * post_shock + debt06 * post_shock +
    post_shock * (white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s +
                    rural_nonfarm_s + poplog_s + under40_s + work_female_s),
  data = g.nn
)
# Create table output
## Estimate variation in outcome 
y0 <- g %>%
  group_by(state) %>%
  summarize(RepVotesMajorPercent=mean(RepVotesMajorPercent,na.rm=TRUE)) %>%
  ungroup() %>%
  summarize(RepVotesMajorPercent=mean(RepVotesMajorPercent,na.rm=TRUE))
y0 <- y0$RepVotesMajorPercent
ysd <- g %>%
  group_by(state) %>%
  summarize(RepVotesMajorPercent=sd(RepVotesMajorPercent,na.rm=TRUE)) %>%
  ungroup() %>%
  summarize(RepVotesMajorPercent=mean(RepVotesMajorPercent,na.rm=TRUE))
ysd <- ysd$RepVotesMajorPercent
#create mean for coal union share
y0u <- g %>%
  filter(!is.na(unionshare)) %>%
  group_by(state) %>%
  summarize(RepVotesMajorPercent=mean(RepVotesMajorPercent,na.rm=TRUE)) %>%
  ungroup() %>%
  summarize(RepVotesMajorPercent=mean(RepVotesMajorPercent,na.rm=TRUE))
y0u <- y0u$RepVotesMajorPercent
ysdu <- g %>%
  filter(!is.na(unionshare)) %>%
  group_by(state) %>%
  summarize(RepVotesMajorPercent=sd(RepVotesMajorPercent,na.rm=TRUE)) %>%
  ungroup() %>%
  summarize(RepVotesMajorPercent=mean(RepVotesMajorPercent,na.rm=TRUE))
ysdu <- ysdu$RepVotesMajorPercent
#Save file output
tab_sens_did <- here("output", "tables", "si_tab_D5_ols_sense.txt")
modelsummary(
  list(m.ols, m.ols.union, m.ols.gfc, m.ols.all),
  vcov = ~ county_state,
  coef_map = coefmap,
  add_rows = data.frame(
    c("Outcome Mean", "Outcome SD", "State Fixed Effects", "Election Fixed Effects"),
    c(round(y0,2), round(ysd,2), "Yes", "Yes"),
    c(round(y0u,2), round(ysdu,2), "Yes", "Yes"),
    c(round(y0,2), round(ysd,2), "Yes", "Yes"),
    c(round(y0u,2), round(ysdu,2), "Yes", "Yes")
  ),
  stars = c("*"=.1,"**"=.05,"***"=.01),
  gof_omit = "IC|RMSE|R2$|Log|Std",
  gof_map = gof_map,
  col.names = c("", "(1)", "(2)", "(3)", "(4)"),
  output = "latex",
  escape = FALSE,
  fmt = 2
) %>%
  group_rows("Treatment and Moderator:", 1, 6) %>%
  group_rows("Time-Varying Covariates:", 7, 12) %>%
  group_rows("Pretreatment Covariates:", 13, 56) %>%
  cat(., file = tab_sens_did)
fix_txt(tab_sens_did)

### Table D.6 ----
#Sensitivity Analysis to Unobserved Confounding
# Regress the treatment on the covariates
D_on_Z <- lm(treated ~ post_shock + coal_bin + post_shock + factor(year) + oilgas_prop + state +
               post_shock * (white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s + rural_nonfarm_s + poplog_s + under40_s + work_female_s),
             g.nn)
D_on_Z.cl <- coeftest(D_on_Z, vcov. = vcovCL(D_on_Z, ~ county_state))
r2dxj.x <- partial_r2(t_statistic = D_on_Z.cl["post_shock:white_s",3], dof = D_on_Z$df.residual)
#regression of the outcome on the treatment and the covariates
r2yxj.dx <- partial_r2(t_statistic = m.ols.cl["post_shock:white_s",3], dof = m.ols$df.residual)
#conduct sensitivity analysis
ols.sens <- sensemakr(estimate = m.ols.cl["treated",][[1]],
                      se = m.ols.cl["treated",][[2]],
                      dof = m.ols$df.residual,
                      r2dxj.x = r2dxj.x,
                      r2yxj.dx = r2yxj.dx,
                      kd = c(5, 10, 15, 20))
ols.sens$bounds$bound_label <- paste0(c("5", "10", "15", "20"), "x White")
summary(ols.sens)
#save table output
tabsensmain <- here("output", "tables", "si_tab_D6_sense.txt")
ols.sens$sensitivity_stats$treatment <- "Coal $\\times$ Post-Shale Shock"
cat(ovb_minimal_reporting(ols.sens), file=tabsensmain)
sense.tab <- readLines(tabsensmain, warn = FALSE)
sense.tab <- sense.tab[-c(1,2,12)]
sense.tab[2] <-  "\\multicolumn{7}{c}{Outcome: Two-Party Republican Presidential Vote Share} \\\\"
sense.tab[3] <- "\\hline"
sense.tab <- c(sense.tab[1], "\\toprule", sense.tab[2:8], "\\bottomrule", sense.tab[9])
cat(sense.tab, file = tabsensmain)
### Fig. D4 ----
dev.off()
png(here("output", "figures", "si_fig_D4_sens.png"), width = 6.5, height = 3, units = "in", res = 300, pointsize = 7)
par(mfrow=c(1,2))
plot(ols.sens, sensitivity.of = "t-value", main = "Sensitivity of t-value", lim=.1, lim.y=.1, label.bump.x	= 0.0075, label.bump.y = 0)
plot(ols.sens, main = "Sensitivity of ATT", lim=.1, lim.y=.1, label.bump.x	= 0.0075, label.bump.y = 0)
dev.off()

## Two-Period Estimates ----
# Estimate with 1972 baseline
m.72to04 <- list()
m.72to04[["2008"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2008), g.nn, subset = (year <= 2004 | year == 2008))
m.72to04[["2012"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2012), g.nn, subset = (year <= 2004 | year == 2012))
m.72to04[["2016"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2016), g.nn, subset = (year <= 2004 | year == 2016))
m.72to04[["2020"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2020), g.nn, subset = (year <= 2004 | year == 2020))
m.72to04[["Pooled"]] <- lm(RepVotesMajorPercent ~ treat * I(year>=2008), g.nn, subset = (year >= 1972))
# Estimate with 1988 baseline
m.88to04 <- list()
m.88to04[["2008"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2008), g.nn, subset = (year %in% c(1988,1992,1996,2000,2004,2008)))
m.88to04[["2012"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2012), g.nn, subset = (year %in% c(1988,1992,1996,2000,2004,2012)))
m.88to04[["2016"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2016), g.nn, subset = (year %in% c(1988,1992,1996,2000,2004,2016)))
m.88to04[["2020"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2020), g.nn, subset = (year %in% c(1988,1992,1996,2000, 2004,2020)))
m.88to04[["Pooled"]] <- lm(RepVotesMajorPercent ~ treat * I(year>=2008), g.nn, subset = (year >= 1988))
#Estimate with 1992 baseline
m.98to04 <- list()
m.98to04[["2008"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2008), g.nn, subset = (year %in% c(1992,1996,2000,2004,2008)))
m.98to04[["2012"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2012), g.nn, subset = (year %in% c(1992,1996,2000,2004,2012)))
m.98to04[["2016"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2016), g.nn, subset = (year %in% c(1992,1996,2000,2004,2016)))
m.98to04[["2020"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2020), g.nn, subset = (year %in% c(1992,1996,2000, 2004,2020)))
m.98to04[["Pooled"]] <- lm(RepVotesMajorPercent ~ treat * I(year>=2008), g.nn, subset = (year >= 1992))
# Estimate with 1996 baseline
m.96to04 <- list()
m.96to04[["2008"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2008), g.nn, subset = (year %in% c(1996,2000,2004,2008)))
m.96to04[["2012"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2012), g.nn, subset = (year %in% c(1996,2000,2004,2012)))
m.96to04[["2016"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2016), g.nn, subset = (year %in% c(1996,2000,2004,2016)))
m.96to04[["2020"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2020), g.nn, subset = (year %in% c(1996,2000,2004,2020)))
m.96to04[["Pooled"]] <- lm(RepVotesMajorPercent ~ treat * I(year>=2008), g.nn, subset = (year >= 1996))
#Estimate with 2000 baseline
m.00to04 <- list()
m.00to04[["2008"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2008), g.nn, subset = (year %in% c(2000,2004,2008)))
m.00to04[["2012"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2012), g.nn, subset = (year %in% c(2000,2004,2012)))
m.00to04[["2016"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2016), g.nn, subset = (year %in% c(2000,2004,2016)))
m.00to04[["2020"]] <- lm(RepVotesMajorPercent ~ treat * I(year==2020), g.nn, subset = (year %in% c(2000,2004,2020)))
m.00to04[["Pooled"]] <- lm(RepVotesMajorPercent ~ treat * I(year>=2008), g.nn, subset = (year >= 2000))
#Create panel
panels <- list(
  "Baseline: 1972-2004" = m.72to04,
  "Baseline: 1988-2004" = m.88to04,
  "Baseline: 1992-2004" = m.98to04,
  "Baseline: 1996-2004" = m.96to04,
  "Baseline: 2000-2004" = m.00to04
)
### Table I2 ----
#Set up table formatting
gm <- gof_map
gm$omit <- TRUE
gm$omit[gm$raw=="nobs"] <- FALSE
#Output table
tab_2x2 <- here("output", "tables", "si_tab_I2_did2x2.txt")
modelsummary(
  panels,
  shape = "rbind",
  vcov = "HC2",
  fmt = 2,
  stars = c("*"=.1,"**"=.05,"***"=.01),
  coef_map = c(
    "treat:I(year == 2008)TRUE"="ATT",
    "treat:I(year == 2012)TRUE"="ATT",
    "treat:I(year == 2016)TRUE"="ATT",
    "treat:I(year == 2020)TRUE"="ATT",
    "treat:I(year == 2004)TRUE"="ATT",
    "treat:I(year >= 2008)TRUE"="ATT"
  ),
  gof_map = gm,
  output = "latex",
  escape = FALSE
) %>%
  add_header_above(c(" " = 1, "Post-Treatment Period" = 5)) %>%
  cat(., file = tab_2x2)
fix_txt(tab_2x2)

## Time-Varying Shale Shocks ----
### Fig. E.4 ----
#State-Specific Shale Shocks Treatment Onset
# Create treatment onset plot
panelview(RepVotesMajorPercent ~ treated_var, data = g.nn, index = c("county_state","year"), 
          axis.lab = "time", xlab = "Election", ylab = "County", 
          background = "white", main = "") %>%
  ggsave(
    .,
    filename = here("output", "figures", "si_fig_E4_state_onset.png"),
    dpi = 300,
    scale = 1.5,
    width = 5.5,
    height = 3
  )

### Table E2 ----
#estimate model
m.var <- fect(RepVotesMajorPercent ~ treated_var + oilgas_prop,
              data = g.nn,
              index = c("county_state", "year"),
              force = "two-way",
              method = "fe",
              se = TRUE,
              nboots = 2000,
              parallel = TRUE,
              na.rm = TRUE)
m.var
#placebo test
m.var.placebo <- fect(RepVotesMajorPercent ~ treated_var + oilgas_prop,
                      data = g.nn,
                      index = c("county_state", "year"),
                      force = "two-way",
                      method = "fe",
                      placeboTest = TRUE,
                      lambda = m.var$lambda.cv, 
                      CV = 0, 
                      se = TRUE,
                      placebo.period = c(-2,0),
                      nboots = 2000,
                      parallel = TRUE,
                      na.rm = TRUE)
m.var.placebo
#Create table output
create_fect_tab(m.var, m.var.placebo, "/si_tab_E2_fect_state.txt")

### Fig. E.5----
var.df <- plot(m.var)$data
p.var <- var.df %>%
  ggplot(aes(x=time,y=ATT,ymin=CI.lower,ymax=CI.upper)) +
  geom_hline(yintercept = 0, color = "grey") +
  geom_vline(xintercept = .5, color = "grey") +
  geom_errorbar(width=0,size=1.25) +
  geom_point(size=3) +
  theme_bw(base_size=14) +
  scale_y_continuous(limits = c(-5, 15)) +
  scale_x_continuous(breaks = seq(-8, 4, 1)) +
  labs(x="Time since the State-Specific Shale Shock Began", y = "Effect on Republican Vote Share (pp)",
       title = "Estimated ATT") +
  theme(plot.title = element_text(hjust=.5),
        panel.grid = element_blank())
#get equiv plot data
var.effectplot <- plot(m.var, type = "equiv")$data
var.minequiv <- var.effectplot$bound[1]
var.equivrange <- var.effectplot$bound[24]
p.varequiv <- var.df %>%
  filter(time <= 0) %>%
  ggplot(aes(x=time,y=ATT,ymin=CI.lower.90,ymax=CI.upper.90)) +
  geom_hline(yintercept=var.minequiv,color="black", lty="dashed") +
  geom_hline(yintercept=-var.minequiv,color="black", lty="dashed") +
  geom_hline(yintercept=var.equivrange,color="red",lty="dashed") +
  geom_hline(yintercept=-var.equivrange,color="red", lty="dashed") +
  geom_hline(yintercept = 0, color = "grey") +
  geom_vline(xintercept = .5, color = "grey") +
  geom_errorbar(width=0,size=1.25) +
  geom_point(size=3) +
  theme_bw(base_size=14) +
  labs(x="Time since the State-Specific Shale Shock Began", y = "Effect on Republican Vote Share (pp)",
       title = "Equivalence Test", lty="",color="") +
  theme(plot.title = element_text(hjust=.5)) +
  scale_y_continuous(limits = c(-5,15)) +
  scale_x_continuous(breaks=seq(-8,0)) +
  annotate("text", label ="p<0.001",x=-4,y=5) +
  scale_color_manual(values=c("red","black")) +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.8, 0.15),
    legend.background = element_rect(fill='transparent'),
    legend.box.background = element_rect(fill='transparent',color="transparent")
  )
var.out <- grid.arrange(p.var, p.varequiv, ncol = 2)
#Create figure output
ggsave(
  var.out,
  filename = here("output", "figures", "si_fig_E5_fect_state.png"),
  width = 6.5,
  height = 2.8,
  dpi = 300,
  scale = 1.5
)

## Coal-to-Gas Ratio Models ----
#estimate model with matched data
price.match <-
  plm(
    formula = RepVotesMajorPercent ~ coal_dd * coalgas_ratio_z + oilgas_prop,
    index = c("county_state","year"),
    data = g.nn
  )
#Cluster standard errors
price.match.cl <- coeftest(price.match, vcov. = plm::vcovHC(price.match, type = "HC2", cluster = "group"))
##unmatched
price.all <-
  plm(
    formula = RepVotesMajorPercent ~ coal_dd * coalgas_ratio_z + oilgas_prop,
    index = c("county_state","year"),
    data = g
  )
#Cluster standard errors
price.all.cl <- coeftest(price.all, vcov. = plm::vcovHC(price.all, type = "HC2", cluster = "group"))
##cross-sectional variance
price.bin <- lm(
  formula = RepVotesMajorPercent ~ coal_dd * coalgas_ratio_z + factor(year) + oilgas_prop + state +
    post_shock * (white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s +
                    rural_nonfarm_s + poplog_s + under40_s + work_female_s),
  data = g
)
#Cluster standard errors
price.bin.cl <- coeftest(price.bin, vcov. = vcovCL(price.bin, type = "HC2", ~ county_state))
#cross-sectional variance plus differential employment
price.emp <- lm(
  formula = RepVotesMajorPercent ~ coal_prop_00_s * coalgas_ratio_z + factor(year) + oilgas_prop + state +
    post_shock * (white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s +
                    rural_nonfarm_s + poplog_s + under40_s + work_female_s),
  data = g
)
#Cluster standard errors
price.emp.cl <- coeftest(price.emp, vcov. = vcovCL(price.emp, type = "HC2", ~ county_state))
##create outcome mean and variance in the treatment and control groups
y0.c <- g.nn %>%
  group_by(county_state) %>%
  summarize(y0 = mean(RepVotesMajorPercent,na.rm=TRUE)) %>%
  ungroup() %>%
  summarize(y0 = mean(y0))
ysd.c <-
  g.nn %>%
  group_by(county_state) %>%
  summarize(y0 = sd(RepVotesMajorPercent,na.rm=TRUE)) %>%
  ungroup() %>%
  summarize(y0 = mean(y0))
y0.c.n <- g %>%
  group_by(county_state) %>%
  summarize(y0 = mean(RepVotesMajorPercent,na.rm=TRUE)) %>%
  ungroup() %>%
  summarize(y0 = mean(y0))
ysd.c.n <-
  g %>%
  group_by(county_state) %>%
  summarize(y0 = sd(RepVotesMajorPercent,na.rm=TRUE)) %>%
  ungroup() %>%
  summarize(y0 = mean(y0))
y0.s <- g %>%
  group_by(state) %>%
  summarize(y0 = mean(RepVotesMajorPercent,na.rm=TRUE)) %>%
  ungroup() %>%
  summarize(y0 = mean(y0))
ysd.s <-
  g %>%
  group_by(state) %>%
  summarize(y0 = sd(RepVotesMajorPercent,na.rm=TRUE)) %>%
  ungroup() %>%
  summarize(y0 = mean(y0))
### Table E3 ----
tab_timevar <- here("output", "tables", "si_tab_E3_timevaryshock.txt")
modelsummary(
  list(price.match, price.all, price.bin, price.emp),
  coef_map = c(
    "(Intercept)" = "Intercept",
    "coal_dd"="Coal",
    "coal_prop_00_s" = "Pre-Shale Shock Coal Employment",
    "post_shock"="Post-Shale Shock",
    "coalgas_ratio_z"="Coal-to-Gas Price Ratio",
    "coal_dd:coalgas_ratio_z"="Coal $\\times$ Post-Shale Shock",
    "coal_prop_00_s:coalgas_ratio_z"="Coal $\\times$ Post-Shale Shock",
    "oilgas_prop"="Hydraulic Fracturing Employment",
    "oilgas_prop_s"="Hydraulic Fracturing Employment",
    "white_s"="White",
    "post_shock:white_s"="White $\\times$ Post-Shale Shock",
    "hispanic_s"="Hispanic",
    "post_shock:hispanic_s"="Hispanic $\\times$ Post-Shale Shock",
    "foreign_s"="Foreign-born",
    "post_shock:foreign_s"="Foreign-born $\\times$ Post-Shale Shock",
    "college_s"="College",
    "post_shock:college_s"="College $\\times$ Post-Shale Shock",
    "Income99pclog_s"="Income per capita (log)",
    "post_shock:income99pclog_s"="Income per capita (log) $\\times$ Post-Shale Shock",
    "poverty_s"="Poverty",
    "post_shock:poverty_s"="Poverty $\\times$ Post-Shale Shock",
    "rural_nonfarm_s"="Rural",
    "post_shock:rural_nonfarm_s"="Rural $\\times$ Post-Shale Shock",
    "poplog_s"="Population (log)",
    "post_shock:poplog_s"="Population (log) $\\times$ Post-Shale Shock",
    "under40_s"="Under 40 years",
    "post_shock:under40_s"="Under 40 years $\\times$ Post-Shale Shock",
    "work_female_s"="Female workforce",
    "post_shock:work_female_s"="Female workforce $\\times$ Post-Shale Shock"
  ),
  vcov = list(
    function(x) plm::vcovHC(x, type="HC2", cluster="group"),
    function(x) plm::vcovHC(x, type="HC2", cluster="group"),
    function(x) sandwich::vcovCL(x, ~ county_state, type = "HC2"),
    function(x) sandwich::vcovCL(x, ~ county_state, type = "HC2")
  ),
  gof_map = gof_map,
  gof_omit = "IC|RMSE|R2$|Log|Std",
  stars = c("*"=.1,"**"=.05,"***"=.01),
  add_rows = data.frame(
    c("Outcome Mean", "Outcome SD",
      "Sample",
      "County Fixed Effects",
      "State Fixed Effects",
      "Election Fixed Effects"),
    c(round(c(y0.c$y0,ysd.c$y0),3), "Matched", "Yes", "No", "Yes"),
    c(round(c(y0.c.n$y0,ysd.c.n$y0),3), "Nation", "Yes", "No", "Yes"),
    c(round(c(y0.s$y0,ysd.s$y0),3), "Nation", "No", "Yes", "Yes"),
    c(round(c(y0.s$y0,ysd.s$y0),3), "Nation", "No", "Yes", "Yes")
  ),
  output = "latex",
  escape = FALSE,
  fmt = 2
) %>%
  group_rows("Treatment and Moderator:", 3,10) %>%
  group_rows("Time-Varying Covariates:", 11,12) %>%
  group_rows("Pretreatment Covariates:", 13,49) %>%
  cat(., file = tab_timevar)
fix_txt(tab_timevar)

## Employment Loss as Treatment ----
### Create matched data using employed data
# Estimate propensity scores
table(g$emptreat,g$adjgroup)
g2000s <- subset(g, year==2000)
fit.emp <- CBPS(update(f.match, emptreat ~ .), data = g2000s, ATT = TRUE)
summary(fit.emp)
out.emp <- matchit(emptreat ~ 1, data = g2000s, distance = fitted(fit.emp), replace = TRUE) 
summary(out.emp)
### Table F1: Balance ----
match.varnames <- c(
  "White",
  "Hispanic",
  "Foreign-born",
  "College",
  "Income per capita (log)",
  "Poverty",
  "Rural",
  "Population (log)",
  "Under 40 years",
  "Female workforce"
)
emp.bal <- rbind(balance(fit.emp)$original, balance(fit.emp)$balanced)
rownames(emp.bal) <- rep(match.varnames, 2)
colnames(emp.bal) <- c("Control", "Treated", "Control", "Treated")
kbl(emp.bal, booktabs = TRUE, digits = 2, escape = FALSE, linesep = "", format="latex") %>%
  add_header_above(c(" "=1,"Mean"=2,"Standardized Mean"=2)) %>%
  group_rows("Original:", 1, 10) %>%
  group_rows("Balanced:", 11, 20) %>%
  cat(., file = here("output", "tables", "si_tab_F1_balance_emp.txt"))
# Prepare data for analysis
df.emp <- match.data(out.emp)
dim(df.emp)
g.emp <- subset(g, fips %in% df.emp$fips)
dim(g.emp)
### Table F2: Balance After Matching----
subset(df.emp, select = c(emptreat, white, hispanic, foreign, college, income99pclog, poverty, rural_nonfarm, poplog, under40, work_female)) %>%
  mutate(treat=ifelse(emptreat==1,"Treatment","Control")) %>%
  dplyr::select(-emptreat) %>%
  dplyr::rename(
    White = white,
    Hispanic = hispanic,
    `Foreign-born` = foreign,
    College = college,
    `Income per capita (log)` = income99pclog,
    Poverty = poverty,
    Rural = rural_nonfarm,
    `Population (log)` = poplog,
    `Under 40 years` = under40,
    `Female workforce` = work_female
  ) %>%
  modelsummary::datasummary_balance(
    formula = White + Hispanic + `Foreign-born` + College + `Income per capita (log)` + Poverty + Rural + `Population (log)` + 
      `Under 40 years` + `Female workforce` ~ treat,
    data = .,
    dinm_statistic = "p.value",
    digits = 2,
    fmt = 2,
    output = "data.frame"
  ) %>%
  kbl(booktabs = TRUE,
      col.names = c("", "Mean", "SD", "Mean", "SD", "Diff.", "$p$"),
      digits = 2,
      format = "latex",
      escape=FALSE,
      linesep="") %>%
  add_header_above(c(" "=1,"Control\n(N=129)"=2,"Treatment\n(N=149)"=2," "=2), escape=TRUE) %>%
  cat(., file = here("output", "tables", "si_tab_F2_matching_emp.txt"))

### Figure F1: Pre-Trends----
emp.pre <- g.emp %>%
  mutate(treat = ifelse(emptreat==1,"Treatment","Control"),
         treat = factor(treat,ordered=TRUE,levels=c("Treatment","Control"))) %>%
  group_by(year, treat) %>%
  summarise(vote = mean(RepVotesMajorPercent, na.rm = TRUE)) %>%
  ggplot(aes(x=year,y=vote,color=treat)) +
  geom_vline(xintercept=2008,color="red") +
  geom_line(aes(lty=treat, color=treat), size = 1.25)+
  geom_point(aes(shape=treat,color=treat), size = 3) +
  labs(color="",shape="",lty="",x = "Election", y = "Two-Party Republican Vote Share") +
  scale_color_grey() +
  theme_bw(base_size = 14) +
  scale_x_continuous(breaks = seq(1968, 2020, 4)) +
  scale_y_continuous(limits = c(30, 80), labels = scales::percent_format(scale=1)) +
  theme(
    panel.grid = element_blank(),
    legend.position = c(.07, .13)
  )
ggsave(
  emp.pre,
  file = here("output", "figures", "si_fig_F1_empmatch_pretrends.png"),
  width = 6.5,
  height = 3, 
  dpi = 300,
  scale = 1.5
)
### Figure F2: Treatment History----
panelview(RepVotesMajorPercent ~ coallayoff_treat, data = g.emp,
          index = c("county_state","year"), 
          axis.lab = "time", xlab = "Election", ylab = "County", 
          background = "white", main = "",
          display.all = TRUE) %>%
  ggsave(
    .,
    filename = here("output", "figures", "si_fig_F2_emp_onset.png"),
    dpi = 300,
    scale = 1.5,
    width = 5.5,
    height = 3
  )
### Figure F3: Dynamic Estimates and Equivalence Test----
fe.emp <- fect(RepVotesMajorPercent ~ coallayoff_treat + oilgas_prop,
               data = g.emp,
               index = c("county_state", "year"),
               force = "two-way",
               method = "fe",
               se = TRUE,
               nboots = 2000,
               parallel = TRUE,
               na.rm = TRUE)
fe.emp

fe.emp.placebo <- fect(RepVotesMajorPercent ~ coallayoff_treat + oilgas_prop,
                       data = g.emp,
                       index = c("fips", "year"),
                       force = "two-way",
                       method = "fe",
                       placeboTest = TRUE,
                       se = TRUE,
                       nboots = 2000,
                       placebo.period = c(-2,0),
                       parallel = TRUE,
                       na.rm = TRUE)
fe.emp.placebo

# Create plot
fe.emp$test.out$tost.equiv.p
emp.out <- plot_fect(fect.model = fe.emp,
                     xmin = -10,
                     equiv.p = "p<0.02",
                     equiv.p.x = -5,
                     xtitle = "Time Since Post-Shale Shock Coal Layoffs")
ggsave(
  emp.out,
  filename = here("output", "figures", "si_fig_F3_fect_emp.png"),
  width = 6.5,
  height = 2.8,
  dpi = 300,
  scale = 1.5
)
### Table F3 -----
create_fect_tab(fe.emp, fe.emp.placebo, "/si_tab_F3_fect_emp.txt")

## Shift-Share Model ------------------------------------------------------
#specify model
f.ss <- RepVotesMajorPercent ~ coal_prop_00_s * nat_coal_z + 
  factor(state) + factor(year) +
  white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s + rural_nonfarm_s + poplog_s + under40_s + work_female_s
#estimate base model
ss <- list()
ss[["(1)"]] <- lm(f.ss, g)
#estimate model with unions
ss[["(2)"]] <- lm(update(f.ss, ~ . + post_shock * scale(unionshare)), g)
#estimate model with fracking controls
ss[["(3)"]] <- lm(update(f.ss, ~ . + post_shock * oilgas_prop), g)
#estimate model with all
ss[["(4)"]] <- lm(update(f.ss, ~ . + post_shock * oilgas_prop + post_shock * scale(unionshare)), g)
#create balancing weights
ss.df <- subset(g, select=c(fips, county_state, coal_prop_00_s, coal_prop_00_s, state, treat, white, hispanic, foreign, college, income99pclog, poverty, rural_nonfarm, poplog, under40, work_female),
                year == 2000 & !is.na(coal_prop_00_s))
ss.df <- unique(ss.df)
dim(ss.df)
fit.ss <- CBPS(update(f.match, coal_prop_00_s ~ .), data = ss.df, ATT = TRUE)
love.plot(fit.ss)
ss.df$weights_ss <- fit.ss$weights
ss.df <- subset(ss.df, select = c(county_state, weights_ss))
g <- left_join(g, ss.df, by = "county_state")
#estimate model with balancing weights
#estimate model without weights
ss[["(5)"]] <- lm(RepVotesMajorPercent ~ coal_prop_00_s * nat_coal_z + post_shock * oilgas_prop + post_shock * scale(unionshare) +
                    factor(year) + factor(state),
                  g)
#add weights
ss[["(6)"]] <- lm(RepVotesMajorPercent ~ coal_prop_00_s * nat_coal_z + post_shock * oilgas_prop + post_shock * scale(unionshare) +
                    factor(year) + factor(state),
                  g, weights = weights_ss)

### Table G1 ----
#create outcome means
ss.outcome.s <-
  g %>%
  filter(!is.na(coal_prop_00_s)) %>%
  group_by(state) %>%
  dplyr::summarize(
    y_avg = mean(RepVotesMajorPercent,na.rm=TRUE),
    y_sd = sd(RepVotesMajorPercent,na.rm=TRUE)
  ) %>%
  ungroup() %>%
  dplyr::summarize(
    y_avg = round(mean(y_avg),3),
    y_sd = round(mean(y_sd),3)
  )
tab_ss <- here("output", "tables", "si_tab_G1_shiftshare.txt")
modelsummary(
  ss,
  stars = c("*"=.1,"**"=.05,"***"=.001),
  vcov = ~ county_state,
  coef_map = c(
    "coal_prop_00_s"="Pre-Shock Coal Employment",
    "nat_coal_z"="National Coal Layoffs",
    "coal_prop_00_s:nat_coal_z"="National Coal Layoffs $\\times$ Post-Shale Shock",
    "oilgas_prop_s"="Hydraulic Fracturing Employment",
    "oilgas_prop"="Hydraulic Fracturing Employment",
    "post_shock:oilgas_prop"="Hydraulic Fracturing Employment $\\times$ Post-Shale Shock",
    "power_prop"="Power Employment",
    "power_prop_s"="Power Employment",
    "power_prop_c"="Power Employment",
    "scale(unionshare)"="Coal Union Share",
    "post_shock:scale(unionshare)"="Coal Union Share $\\times$ Post-Shale Shock",
    "white_s"="White",
    "post_shock:white_s"="White $\\times$ Post-Shale Shock",
    "hispanic_s"="Hispanic",
    "post_shock:hispanic_s"="Hispanic $\\times$ Post-Shale Shock",
    "foreign_s"="Foreign-born",
    "post_shock:foreign_s"="Foreign-born $\\times$ Post-Shale Shock",
    "college_s"="College",
    "post_shock:college_s"="College $\\times$ Post-Shale Shock",
    "income99pclog_s"="Income per capita (log)",
    "post_shock:income99pclog_s"="Income per capita (log) $\\times$ Post-Shale Shock",
    "poverty_s"="Poverty",
    "post_shock:poverty_s"="Poverty $\\times$ Post-Shale Shock",
    "rural_nonfarm_s"="Rural",
    "post_shock:rural_nonfarm_s"="Rural $\\times$ Post-Shale Shock",
    "poplog_s"="Population (log)",
    "post_shock:poplog_s"="Population (log) $\\times$ Post-Shale Shock",
    "under40_s"="Under 40 years",
    "post_shock:under40_s"="Under 40 years $\\times$ Post-Shale Shock",
    "work_female_s"="Female workforce",
    "post_shock:work_female_s"="Female workforce $\\times$ Post-Shale Shock"
  ),
  add_rows = data.frame(
    c("Outcome Mean", "Outcome SD", "State Fixed Effects",
      "Election Fixed Effects", "Covariate Balancing Propensity Score"),
    c(ss.outcome.s$y_avg, ss.outcome.s$y_sd, "Yes", "Yes", "No"),
    c(ss.outcome.s$y_avg, ss.outcome.s$y_sd, "Yes", "Yes", "No"),
    c(ss.outcome.s$y_avg, ss.outcome.s$y_sd, "Yes", "Yes", "No"),
    c(ss.outcome.s$y_avg, ss.outcome.s$y_sd, "Yes", "Yes", "No"),
    c(ss.outcome.s$y_avg, ss.outcome.s$y_sd, "Yes", "Yes", "No"),
    c(ss.outcome.s$y_avg, ss.outcome.s$y_sd, "No", "Yes", "Yes")
  ),
  gof_omit = "IC|RMSE|R2$|Log|Std",
  gof = gof_map,
  output = "latex",
  escape = FALSE,
  fmt = 2
) %>%
  group_rows("Treatment and Moderator:", 1, 6) %>%
  group_rows("Time-Varying Covariates:", 7, 14) %>%
  group_rows("Pretreatment Covariates:", 15, 33) %>%
  cat(., file = tab_ss)
fix_txt(tab_ss)
###Figure G1-----
#Estimate dynamic model
ss.ct <-
  lm(RepVotesMajorPercent ~ coal_prop_00_s * factor(year) +
       white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s + rural_nonfarm_s + poplog_s + under40_s + work_female_s +
       oilgas_prop + factor(state),
     g) %>%
  lmtest::coeftest(., vcov. = sandwich::vcovCL(., ~ county_state))
ss.ct <- broom::tidy(ss.ct, conf.int = TRUE)
ss.ct <- subset(ss.ct, grepl("coal_prop_00_s:factor", term))
ss.ct$term <- gsub("coal_prop_00_s:factor\\(year\\)", "", ss.ct$term)
ss.ct$term <- as.numeric(ss.ct$term)
preshockavg <- mean(subset(ss.ct, term < 2008)$estimate)
ss.ct <-
  rbind(
    ss.ct,
    data.frame(
      term = 1976,
      estimate = 0,
      std.error = NA,
      statistic = NA,
      p.value = NA,
      conf.low = NA,
      conf.high = NA
    )
  )
p.ss.ct <-
  ss.ct %>%
  ggplot(aes(x=term,y=estimate,ymin=conf.low,ymax=conf.high)) +
  geom_vline(xintercept = 2008, color = "grey") +
  geom_hline(yintercept = 0, color = "grey") +
  geom_errorbar(width = 0, size = 1.25) +
  geom_point(size = 3) +
  geom_hline(yintercept = preshockavg, color = "grey", lty = "dashed") +
  annotate("text", x = 2016, y = preshockavg - .25, label = "Pre-Shock Average") +
  theme_bw(base_size = 14) +
  scale_x_continuous(breaks = seq(1976, 2020, 4)) +
  labs(x="Election", y="Effect on Republican Vote Share") +
  theme(
    panel.grid = element_blank()
  )
p.ss.ct
ggsave(
  p.ss.ct,
  filename = here("output", "figures", "si_fig_G1_ss_trends.png"),
  width = 6.5,
  height = 3,
  scale = 1.5,
  dpi = 300
)
#Heterogeneous Treatment Effects ----
## Unions ----
#Estimate model controlling for unionization levels
m.union <- fect(RepVotesMajorPercent ~ treated_var + oilgas_prop + unionshare,
                data = g.nn,
                index = c("county_state", "year"),
                force = "two-way",
                method = "fe",
                CV = TRUE,
                se = TRUE,
                nboots = 2000,
                parallel = TRUE,
                na.rm = TRUE)
m.union
m.union.placebo <- fect(RepVotesMajorPercent ~ treated_var + oilgas_prop + unionshare,
                        data = g.nn,
                        index = c("county_state", "year"),
                        force = "two-way",
                        method = "fe",
                        placeboTest = TRUE,
                        se = TRUE,
                        placebo.period = c(-3, 0),
                        nboots = 2000,
                        parallel = TRUE,
                        na.rm = TRUE)
### Table H1----
att.tab <-
  rbind(
    m.union$est.avg,
    m.union$est.avg.unit,
    m.union$est.beta,
    m.union.placebo$est.placebo[,c(1:5)]
  )
colnames(att.tab) <- c("Estimate", "S.E.", "CI$_{2.5\\%}$", "CI$_{97.5\\%}$", "$p$-value")
rownames(att.tab) <- c("Observations equally weighted", "Units equally weighted",
                       "Hydraulic Fracturing Employment", "Unionized Coal Production Share",
                       "-3 to 0 election interval")
kbl(att.tab, booktabs=TRUE, digits = 2, linesep = "",
    escape = FALSE, format="latex") %>%
  group_rows("ATT:", 1, 2) %>%
  group_rows("Covariates:", 3, 4) %>%
  group_rows("Placebo Tests:", 5, 5) %>%
  cat(., file = here("output", "tables", "si_tab_H1_union_fect.txt"))

### Figure H3----
p.union <- plot_fect(fect.model = m.union, equiv.p = "p<0.04", ymax = 15, equiv.p.x=-2)

ggsave(
  p.union,
  filename = here("output", "figures", "si_fig_H3_fect_union.png"),
  width = 6.5,
  height = 2.8,
  dpi = 300,
  scale = 1.5
)

### Heterogeneity ----
#subset the data to pre and post
df.treat <- subset(g.nn, treat == 1)
df.treat <- as.data.frame(df.treat)

#### Figure H5 ----
#check linearity assumption
raw <-
  interflex(
    estimator = "raw",
    data = df.treat,
    Y = "RepVotesMajorPercent",
    D = "post_shock",
    X = "unionshare",
    Z = c("oilgas_prop"),
    FE = c("county_state", "year"),
    na.rm = TRUE,
    treat.type = "discrete",
    vcov.type = "cluster",
    cl = "county_state",
    show.grid = FALSE,
    theme.bw = TRUE
  )
raw <-
  raw +
  labs(x="Proportion of Unionized Coal Production, County-Level", 
       y = "Republican Presidential Vote Share") +
  theme_bw(base_size = 14) +
  theme(panel.grid = element_blank())
raw
ggsave(
  raw,
  filename = here("output", "figures", "si_fig_H5_unionlin.png"),
  width = 5,
  height = 3.5,
  dpi = 300,
  scale = 1.5
)
#binning estimator
bin <-
  interflex(
    estimator = "binning",
    data = df.treat,
    Y = "RepVotesMajorPercent",
    D = "post_shock",
    X = "unionshare",
    Z = c("oilgas_prop"),
    FE = c("county_state", "year"),
    na.rm = TRUE,
    treat.type = "discrete",
    vcov.type = "cluster",
    cl = "county_state",
    theme.bw = TRUE,
    show.grid = FALSE,
    bin.labs = FALSE,
    xlab = "Proportion of Unionized Coal Production, County-Level",
    ylab = "Effect on Republican Vote Share"
  )
bin

#kernel estimator
kern <-
  interflex(
    estimator = "kernel",
    data = df.treat,
    Y = "RepVotesMajorPercent",
    D = "post_shock",
    X = "unionshare",
    Z = c("oilgas_prop"),
    FE = c("county_state", "year"),
    na.rm = TRUE,
    treat.type = "discrete",
    vcov.type = "cluster",
    cl = "county_state",
    theme.bw = TRUE,
    show.grid = FALSE,
    bin.labs = FALSE,
    nboots = 2000,
    xlab = "Proportion of Unionized Coal Production, County-Level",
    ylab = "Effect on Republican Vote Share")
kern

union.p <- gridExtra::grid.arrange(bin$figure, kern$figure,ncol=2)
union.p


## Transitional Assistance / POWER Grants ----
#Create counties that received funds
g.funds.fips <- subset(g.nn, funds_total > 0)
g.funds <- subset(g.nn, fips %in% g.funds.fips$fips)
`%notin%` <- Negate(`%in%`)
g.nofunds <- subset(g.nn, fips %notin% g.funds.fips$fips)
# Set up covariates
g.nn.funds <- g.nn
g.nn.funds <- g.nn.funds %>%
  group_by(fips) %>%
  mutate(across(funds_arc:grants_cum, .fns = list("sd" = ~ scale(.x)[,1]), .names = "{.col}_{.fn}"))
g.nn.funds <- g.nn.funds %>%
  mutate(across(funds_arc:grants_cum, .fns = list("obama" = ~ case_when(year > 2016 ~ 0, year <= 2016 ~ .x)), .names = "{.col}_{.fn}"))
g.nn.funds <- g.nn.funds %>%
  group_by(fips) %>%
  mutate(across(funds_arc_obama:grants_cum_obama, .fns = list("sd" = ~ scale(.x)[,1]), .names = "{.col}_{.fn}"))

# Estimate models
m.power <- list()
m.power[["(1)"]] <- felm(
  RepVotesMajorPercent ~ treated * funds_total_sd | county_state + year | 0 | county_state,
  g.nn.funds
)
m.power[["(2)"]] <- felm(
  RepVotesMajorPercent ~ treated * funds_total_cum_sd | county_state + year | 0 | county_state,
  g.nn.funds
)
m.power[["(3)"]] <- felm(
  RepVotesMajorPercent ~ treated * grants | county_state + year | 0 | county_state,
  subset(g.nn.funds, !is.na(funds_total_sd))
)
m.power[["Yes"]] <- felm(
  RepVotesMajorPercent ~ treated | county_state + year | 0 | county_state,
  data = g.funds
)
m.power[["No"]] <- felm(
  RepVotesMajorPercent ~ treated | county_state + year | 0 | county_state,
  data = g.nofunds
)
### Table H2 ----
tab_power <- here("output", "tables", "si_tab_H2_tab_power.txt")
modelsummary(
  m.power,
  stars = c("*"=.1,"**"=.05,"***"=.01),
  gof_map = gof_map,
  coef_map = c(
    "treated"="Treatment",
    "funds_total_sd" = "POWER Funding",
    "treated:funds_total_sd" = "Treatment $\\times$ POWER Funding",
    "funds_total_cum_sd" = "Cumulative POWER Funding",
    "treated:funds_total_cum_sd" = "Treatment $\\times$ Cumulative POWER Funding",
    "grants"="POWER Grants",
    "treated:grants"="Treatment $\\times$ POWER Grants"
  ),
  gof_omit = "R2|IC|RMSE|St",
  fmt = 2,
  escape = FALSE,
  output = "latex",
  add_rows = data.frame(
    c("County Fixed Effects", "Election Year Fixed Effects"),
    c("Yes", "Yes"),
    c("Yes", "Yes"),
    c("Yes", "Yes"),
    c("Yes", "Yes"),
    c("Yes", "Yes")
  )
) %>%
  add_header_above(c(" " = 4, "Funding" = 2)) %>%
  cat(., file = tab_power)
fix_txt(tab_power)

# Falsification Tests----
## Estimate models
# Iron ore
plm.iron <- 
  plm(
    formula = RepVotesMajorPercent ~ iron_falsify * post_shock,
    index = c("county_state", "year"),
    model = "within",
    data = subset(g, coal_dd == 0)
  )
plm.iron.cl <- coeftest(plm.iron, plm::vcovHC(plm.iron, type = "HC2", cluster = "group"))
plm.iron.cl
# Silver mine
plm.silver <- 
  plm(
    formula = RepVotesMajorPercent ~ silver_falsify * post_shock,
    index = c("county_state", "year"),
    model = "within",
    data = subset(g, coal_dd == 0)
  )
plm.silver.cl <- coeftest(plm.silver, plm::vcovHC(plm.silver, type = "HC2", cluster = "group"))
plm.silver.cl
# Sandstone mines
plm.sandstone <- 
  plm(
    formula = RepVotesMajorPercent ~ sandstone_falsify * post_shock,
    index = c("county_state", "year"),
    model = "within",
    data = subset(g, coal_dd == 0)
  )
plm.sandstone.cl <- coeftest(plm.sandstone, plm::vcovHC(plm.sandstone, type = "HC2", cluster = "group"))
plm.sandstone.cl
# General layoffs
plm.layoff <- 
  plm(
    formula = RepVotesMajorPercent ~ layoff * post_shock,
    index = c("county_state", "year"),
    data = subset(g, coal_dd == 0)
  )
plm.layoffs.cl <- coeftest(plm.layoff, plm::vcovHC(plm.layoff, type = "HC2", cluster = "group"))
plm.layoffs.cl

## Figure J2: Pre-Trends ----
# Create theme for pre-trend plots
false_theme <- list(
  geom_vline(xintercept=2008, color="red"),
  geom_line(aes(x=year,y=RepVotesMajorPercent,color=test,lty=test), size = 1.25),
  geom_point(aes(x=year,y=RepVotesMajorPercent,color=test,shape=test), size = 3),
  theme_bw(base_size = 14),
  scale_color_grey(),
  theme(
    panel.grid = element_blank(),
    legend.position = c(.88, .19),
    legend.margin = margin(0),
    legend.text.align = 0,
    legend.text = element_text(size=8),
    legend.background = element_rect(fill='transparent'),
    legend.box.background = element_rect(fill='transparent',color="transparent")
  ),
  scale_x_continuous(breaks = seq(1972, 2020, 8)),
  scale_y_continuous(limits = c(30, 80), labels = scales::percent_format(scale=1)),
  labs(x = "Year", y = "Two-Party Republican Vote Share", shape = "", lty = "", color = "")
)
# Create pre-trends plot for silver mining
silver.data <-
  g %>%
  dplyr::mutate(
    test = dplyr::case_when(
      silver_falsify == 1 ~ "Silver",
      adjgroup == "Coal" ~ "Coal",
      silver_falsify == 0 & is.na(adjgroup) ~ "Neither"
    )
  ) %>%
  dplyr::filter(!is.na(test)) %>%
  dplyr::group_by(year, test) %>%
  dplyr::summarise(RepVotesMajorPercent = mean(RepVotesMajorPercent,na.rm=TRUE))
p.silver <-
  silver.data %>%
  ggplot() +
  false_theme +
  ggtitle("Silver Mine")
p.silver

# Create pre-trends plot for iron ore
iron.data <-
  g %>%
  dplyr::mutate(
    test = dplyr::case_when(
      iron_falsify == 1 ~ "Iron",
      adjgroup == "Coal" ~ "Coal",
      iron_falsify == 0 & is.na(adjgroup) ~ "Neither"
    )
  ) %>%
  dplyr::filter(!is.na(test)) %>%
  dplyr::group_by(year, test) %>%
  dplyr::summarise(RepVotesMajorPercent = mean(RepVotesMajorPercent, na.rm = TRUE))
p.iron <-
  iron.data %>%
  ggplot() +
  false_theme +
  ggtitle("Iron Mine")
p.iron

# Create pre-trends plot for sandstone mining
sand.data <-
  g %>%
  dplyr::mutate(
    test = dplyr::case_when(
      sandstone_falsify == 1 ~ "Sandstone",
      adjgroup == "Coal" ~ "Coal",
      sandstone_falsify == 0 & is.na(adjgroup) ~ "Neither"
    )
  ) %>%
  dplyr::filter(!is.na(test)) %>%
  dplyr::group_by(year, test) %>%
  dplyr::summarise(RepVotesMajorPercent = mean(RepVotesMajorPercent,na.rm=TRUE))
p.sand <-
  sand.data %>%
  ggplot() +
  false_theme +
  ggtitle("Sandstone Mine")
p.sand

# Create pre-trends plot for general layoffs
layoff.data <-
  g %>%
  dplyr::mutate(
    test = dplyr::case_when(
      layoff == 1 & coal_dd == 0 ~ "Layoffs",
      adjgroup == "Coal" ~ "Coal",
      sandstone_falsify == 0 & is.na(adjgroup) ~ "Neither"
    )
  ) %>%
  dplyr::filter(!is.na(test)) %>%
  dplyr::group_by(year, test) %>%
  dplyr::summarise(RepVotesMajorPercent = mean(RepVotesMajorPercent,na.rm=TRUE))
p.layoff <-
  layoff.data %>%
  ggplot() +
  false_theme +
  ggtitle("Layoffs")
p.layoff

# Combine pre-trend plots
false.pretrends <- grid.arrange(p.iron,p.sand,p.silver,p.layoff,nrow=2,ncol=2)
ggsave(
  false.pretrends,
  filename = here("output", "figures", "si_fig_J2_false_pretrends.png"),
  width = 6.5,
  height = 5,
  scale = 1.5,
  dpi = 300
)

## Figure J3: Anthracite Falsification Test ----
# Subset to coal mine counties
ant <- subset(g, (coal_dd == 1 | ant_treat_emp == 1))
# Define treatment status
ant$treated <- ifelse(ant$year >= ant$break_year & ant$ant_treat_emp == 0, 1, 0)
# Estimate DiD with FEct
m.ant <- fect(RepVotesMajorPercent ~ treated + oilgas_prop,
              data = ant,
              index = c("county_state", "year"),
              force = "two-way",
              method = "fe",
              CV = TRUE,
              se = TRUE,
              nboots = 2000,
              parallel = TRUE,
              na.rm = TRUE)
m.ant
# Create plot
p.ant <- plot_fect(m.ant, equiv.p = "p<0.001", ymax=15, equiv.p.x=-3)
ggsave(
  p.ant,
  filename = here("output", "figures", "si_fig_J3_fect_ant.png"),
  width = 6.5,
  height = 2.8,
  dpi = 300,
  scale = 1.5
)

## Model With Continuous Employment Treatment ----
m.emp <- lm(RepVotesMajorPercent ~ coal_prop_00_c * post_shock + oilgas_prop + factor(year) + factor(county_state), g.nn)
m.emp.cl <- coeftest(m.emp, vcov. = vcovCL(m.emp, ~ county_state))
# Check N
subset(g, year==2000 & coal_prop_00_c > 1) %>% nrow()

## Figure J1: Summary of Falsification Tests ----
false.df <- rbind(m.basic$est.avg[,c(1:2)], m.ant$est.avg[,c(1:2)],
                  fe.emp$est.avg[,c(1:2)],
                  m.emp.cl["coal_prop_00_c:post_shock",c(1:2)],
                  plm.iron.cl[2,c(1,2)], plm.silver.cl[2,c(1,2)], plm.sandstone.cl[2,c(1,2)], plm.layoffs.cl[3,c(1,2)])
false.df <- as.data.frame(false.df)
#determine sample size
m.basic$Nco
m.basic$Ntr
dim(g.nn)
m.ant$Nco
m.ant$Ntr
dim(ant)
fe.emp$Nco
fe.emp$Ntr
dim(g.emp)
n_distinct(g.nn$county_state)
summary(plm.iron)
n_distinct(subset(g, iron_falsify==1)$county_state)
n_distinct(subset(g, iron_falsify==0 & coal_dd == 0)$county_state)
summary(plm.sandstone)
n_distinct(subset(g, sandstone_falsify==1)$county_state)
n_distinct(subset(g, sandstone_falsify==0 & coal_dd == 0)$county_state)
summary(plm.silver)
n_distinct(subset(g, silver_falsify==1)$county_state)
n_distinct(subset(g, silver_falsify==0 & coal_dd == 0)$county_state)
summary(plm.layoff)
n_distinct(subset(g, layoff==1 & coal_dd == 0)$county_state)
layofffips <- subset(g, layoff==1 & coal_dd == 0)$fips %>% unique()
`%notin%` <- Negate(`%in%`)
n_distinct(subset(g, fips %notin% layofffips & coal_dd == 0)$county_state)
# Create model names
modelnames <- c("Coal \nvs.\nMatched",
                "Vulnerable\nvs.\nUnexposed",
                "Coal Layoffs\nvs.\nMatched",
                "+1SD Shock\nExposure",
                "Iron\nvs.\nNation",
                "Sandstone\nvs.\nNation",
                "Silver\nvs.\nNation",
                "+1SD Layoffs\nvs.\nNation")
false.df$resource <- modelnames
names(false.df) <- c("est", "se", "resource")
false.df$lb <- false.df$est + qnorm(0.025) * false.df$se
false.df$ub <- false.df$est + qnorm(0.975) * false.df$se
false.df$test <- c(rep("Coal Country (Treatment)",4), rep("Other Mining Counties (Placebo)", 3),
                   "Counties with Non-Coal Layoffs (Placebo)")
false.df$test <- factor(false.df$test, ordered = TRUE, levels = c("Coal Country (Treatment)",
                                                                  "Other Mining Counties (Placebo)",
                                                                  "Counties with Non-Coal Layoffs (Placebo)"))
false.df$resource <- factor(false.df$resource, ordered = TRUE, levels = modelnames)
# Create plot
false.plot <- false.df %>%
  ggplot() +
  geom_hline(yintercept = 0, color="grey")+
  geom_pointrange(aes(x = resource, y = est, ymin = lb, ymax = ub), fatten = 4.5,
                  linewidth = 1.5) +
  theme_bw(base_size = 14) +
  facet_grid(~test, scales="free_x", space = "free_x", labeller = label_wrap_gen(20)) +
  labs(color ="", shape ="", x="Treatment/Placebo", y = "Effect on Republican Vote Share (pp)", title= "ATT of Shale Shock")+
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.9, 0.15),
    plot.title = element_text(hjust =.5)
  )
false.plot
ggsave(
  false.plot,
  filename = here("output", "figures", "si_fig_J1_falsetest.png"),
  width = 6.5,
  height = 3,
  dpi = 300,
  scale = 1.5
  )


## Table J2 ----
false_models <- list(plm.iron.cl, plm.silver.cl, plm.sandstone.cl, plm.layoffs.cl)
names(false_models) <- c("Iron", "Silver", "Sandstone", "Layoffs")

tab_falsetest <- here("output", "tables", "si_tab_J2_falsificationtest.txt")
modelsummary(
  false_models,
  stars = c("*"=.1,"**"=.05,"***"=.01),
  gof_map = NA,
  add_rows = data.frame(
    c("N", "Adjusted $R^2$", "County Fixed Effects", "Election Fixed Effects"),
    c(format(nobs(plm.iron), nsmall = 0),format(round(summary(plm.iron)$r.squared[[2]],2), nsmall=2), "Yes", "Yes"),
    c(format(nobs(plm.silver), nsmall = 0),format(round(summary(plm.silver)$r.squared[[2]],2), nsmall=2), "Yes", "Yes"),
    c(format(nobs(plm.sandstone), nsmall = 0),format(round(summary(plm.sandstone)$r.squared[[2]],2), nsmall=2), "Yes", "Yes"),
    c(format(nobs(plm.layoff), nsmall = 0),format(round(summary(plm.layoff)$r.squared[[2]],2), nsmall=2),"Yes", "Yes")
  ),
  coef_map = c(
    "post_shock"="Post-Shale Shock",
    "treat:post_shock"="Coal x Post-Shale Shock",
    "coal_dd:post_shock"="Coal x Post-Shale Shock",
    "iron_falsify:post_shock"="Iron Mine x Post-Shale Shock",
    "silver_falsify:post_shock"="Silver Mine x Post-Shale Shock",
    "sandstone_falsify:post_shock"="Sandstone Mine x Post-Shale Shock",
    "layoff:post_shock"="Non-Coal Layoffs x Post-Shale Shock",
    "layoff" = "Non-Coal Layoffs"
  ),
  escape = FALSE,
  output = "latex",
  fmt = 2
) %>%
  cat(., file = tab_falsetest)

false.tab <- readLines(tab_falsetest)
false.tab <- false.tab[-c(1:2, 24, 26)]
false.tab <- c(false.tab[1:16], "\\midrule", false.tab[17:22])
cat(false.tab, file = tab_falsetest)

## Table J1: Covariate Balance----
table(g$iron_falsify,g$layoff)
table(g$silver_falsify,g$layoff)
table(g$sandstone_falsify,g$layoff)
table(g$coal_dd,g$layoff)
layoff_fips<-subset(g,layoff==1&coal_dd==0)
g$layoff_treat <- ifelse(g$fips%in%layoff_fips$fips,1,0)
# Define groups
g <- g %>%
  mutate(
    samplefalse = case_when(
      iron_falsify==1~"Iron",
      coal_dd==1~"Coal",
      silver_falsify==1~"Silver",
      sandstone_falsify==1~"Sandstone",
      layoff_treat==1~"Layoff",
      T ~ NA_character_
    )
  )
#summary stats for demographics
dem <- g %>%
  filter(!is.na(samplefalse)) %>%
  rename(workfemale=work_female,ruralnonfarm=rural_nonfarm) %>%
  group_by(samplefalse) %>%
  summarise(
    across(
      c(white, hispanic, foreign, college, income99pclog, poverty, ruralnonfarm, poplog, under40, debt06, workfemale),
      .fns = list("mean"=~mean(.x,na.rm=TRUE),"sd"=~sd(.x,na.rm=TRUE))
    )
  ) %>%
  pivot_longer(cols = c(white_mean:workfemale_sd)) %>%
  separate(col=name,into=c("var","stat"), sep="_") %>%
  pivot_wider(names_from=c(samplefalse,stat))
#summary stats for time-varying
tvar <- g %>%
  filter(!is.na(samplefalse)) %>%
  rename(gasdist=NewGasDistance_s,oilgasprop=oilgas_prop,coalprop=coal_prop,powerprop=power_prop) %>%
  group_by(samplefalse) %>%
  summarise(
    across(
      c(gasdist,oilgasprop,coalprop,powerprop,RepVotesMajorPercent),
      .fns = list("mean"=~mean(.x,na.rm=TRUE),"sd"=~sd(.x,na.rm=TRUE))
    )
  ) %>%
  pivot_longer(cols = c(gasdist_mean:RepVotesMajorPercent_sd)) %>%
  separate(col=name,into=c("var","stat"), sep="_") %>%
  pivot_wider(names_from=c(samplefalse,stat))
#total number
D_coal <- subset(g,coal_dd==1,select=county_state) %>% n_distinct()
N_coal <- subset(g,coal_dd==1) %>% nrow()
D_iron <- subset(g,samplefalse=="Iron",select=county_state) %>% n_distinct()
N_iron <- subset(g,samplefalse=="Iron") %>% nrow()
D_layoff <- subset(g,samplefalse=="Layoff",select=county_state) %>% n_distinct()
N_layoff <- subset(g,samplefalse=="Layoff") %>% nrow()
D_sand <- subset(g,samplefalse=="Sandstone",select=county_state) %>% n_distinct()
N_sand <- subset(g,samplefalse=="Sandstone") %>% nrow()
D_silver <- subset(g,samplefalse=="Silver",select=county_state) %>% n_distinct()
N_silver <- subset(g,samplefalse=="Silver") %>% nrow()
N_info <-
  data.frame(
    var = c("Units", "Total $N$"),
    Coal_mean = c(D_coal, N_coal),
    Coal_sd = c(NA,NA),
    Iron_mean = c(D_iron,N_iron),
    Iron_sd = c(NA,NA),
    Layoff_mean = c(D_layoff,N_layoff),
    Layoff_sd = c(NA,NA),
    Sandstone_mean = c(D_sand,N_sand),
    Sandstone_sd = c(NA,NA),
    Silver_mean = c(D_silver,N_silver),
    Silver_sd = c(NA,NA)
  )
#rearrange dem to match tvar and n_info
falsesum <-
  bind_rows(
    digitsByRows(dem[,c(2:11)], 2),
    digitsByRows(tvar[,c(2:11)], 2),
    digitsByRows(N_info[,c(2:11)], 0)
  )
#bind rows together
falsesum$var <- c("White", "Hispanic", "Foreign-born", "College", "Income per capita (log)",
                  "Poverty", "Rural", "Population (log)", "Under 40 years", "Debt-to-income",
                  "Female workforce", "Gas plant distance (1992--2020)", "Hydraulic fracturing employment", "Coal employment",
                  "Power employment", "Two-party Republican vote share", "Counties", "$N$")
falsesum <- falsesum[,c(11,1:4,7:10,5:6)]
names(falsesum) <- c("", "Mean", "SD", "Mean", "SD", "Mean", "SD", "Mean", "SD","Mean","SD")
kbl(falsesum, booktabs = TRUE, escape = FALSE, format = "latex", linesep = "", digits=2) %>%
  add_header_above(c(" "=1, "Coal"=2,"Iron"=2, "Sandstone"=2,"Silver"=2, "Layoffs"=2)) %>%
  group_rows("Pretreatment:", 1, 11) %>%
  group_rows("Time-Varying:", 12, 16) %>%
  group_rows("Observations:", 17, 18) %>%
  cat(., file = here("output", "tables", "si_tab_J1_false_balance.txt"))

# Fig. 4: Moderating Effect of Visibility ----
## Estimate CBPS----
g$treatvis <- with(g, treat * NewGasDistance_s)
gvis <- subset(g, year == 2008 & !is.na(treatvis))
fitvis <- CBPS(update(f.match, treatvis ~ .), ATT = TRUE, data = gvis, method = "exact")
p.fitvis.bal <-
  love.plot(fitvis,
            var.names = c("white"="White", "hispanic"="Hispanic", "foreign"="Foreign-Born",
                          "college"="College", "income99pclog"="Income per capita (log)",
                          "poverty"="Poverty", "rural_nonfarm"="Rural",
                          "poplog"="Population (log)", "under40"="Under 40 years",
                          "work_female"="Female workforce"),
            colors = c("darkgrey", "black"),
            shapes = c(19, 17)) +
  theme_bw(base_size = 14) +
  labs(title = "") +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.9, .1)
  )
### Figure K.2----
## Covariate Balance by Continuous Treatment Status in the Adjusted and Un-adjusted Samples
ggsave(
  p.fitvis.bal,
  filename = here("output", "figures", "si_fig_K2_fitvis_balance.png"),
  scale = 1.5,
  dpi = 300,
  width = 6,
  height = 4
)
#add weights to data
gvis$weights_vis <- fitvis$weights
gvis <- subset(gvis, select = c(county_state, weights_vis))
g <- left_join(g, gvis, by = "county_state")

## Estimate Models ----
#cross-sectional variation by state
m.vis.cs <-
  lm(
    formula = RepVotesMajorPercent ~ treat * post_shock * NewGasDistance_s + factor(year) + state +
      oilgas_prop + post_shock * unionshare +
      white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s + rural_nonfarm_s + poplog_s + under40_s + work_female_s,
    data = subset(g, year >= 1992 & !is.na(treat))
  )
m.vis.cs.cl <- coeftest(m.vis.cs, vcov. = vcovCL(m.vis.cs, ~ county_state))
m.vis.cs.cl
#add control for power employment
m.vis.cs.power <-
  lm(
    formula = RepVotesMajorPercent ~ treat * post_shock * NewGasDistance_s + factor(year) + state +
      oilgas_prop + power_prop + post_shock * unionshare +
      white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s + rural_nonfarm_s + poplog_s + under40_s + work_female_s,
    data = subset(g, year >= 1992 & !is.na(treat))
  )
m.vis.cs.power.cl <- coeftest(m.vis.cs.power, vcov. = vcovCL(m.vis.cs.power, ~ county_state))
m.vis.cs.power.cl
#interact moderator with covariates
m.inter <-
  lm(
    formula = RepVotesMajorPercent ~ treat * post_shock * NewGasDistance_s + factor(year) + state +
      NewGasDistance_s * (oilgas_prop + power_prop + post_shock * unionshare +
                            white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s + rural_nonfarm_s + poplog_s + under40_s + work_female_s),
    data = subset(g, year >= 1992 & !is.na(treat))
  )
m.inter.cl <- coeftest(m.inter, vcov. = vcovCL(m.inter, ~ county_state))
m.inter.cl
#dif-in-dif
m.vis.dd <- lm(RepVotesMajorPercent ~ treat * post_shock * NewGasDistance_c + oilgas_prop + power_prop + unionshare * post_shock +
                 factor(county_state) + factor(year),
               data = subset(g, year>=1992 & !is.na(treat)))
m.vis.dd.cl <- coeftest(m.vis.dd, vcov. = vcovCL(m.vis.dd, ~ county_state))
tail(m.vis.dd.cl)
#dif-in-dif with interacted moderator
m.inter.dd <- lm(RepVotesMajorPercent ~ treat * post_shock * NewGasDistance_c + NewGasDistance_c * (oilgas_prop + power_prop + unionshare * post_shock) +
                   factor(county_state) + factor(year),
                 data = subset(g, year>=1992 & !is.na(treat)))
m.inter.dd.cl <- coeftest(m.inter.dd, vcov. = vcovCL(m.inter.dd, ~ county_state))
tail(m.inter.dd.cl)
#weighted dif-in-dif
m.vis.dd.pi <- lm(RepVotesMajorPercent ~ treat * post_shock * NewGasDistance_c + oilgas_prop + power_prop + unionshare * post_shock +
                    factor(county_state) + factor(year),
                  data = g,
                  weights = weights_vis
)
m.vis.dd.pi.cl <- coeftest(m.vis.dd.pi, vcov. = vcovCL(m.vis.dd.pi, ~ county_state))
tail(m.vis.dd.pi.cl)
#weighted cross sectional model
m.vis.cs.wt <-
  lm(
    formula = RepVotesMajorPercent ~ treat * post_shock * NewGasDistance_s + factor(year) + state +
      oilgas_prop + power_prop + post_shock * unionshare,
    data = subset(g, year >= 1992 & !is.na(treat)),
    weights = weights_vis
  )
m.vis.cs.wt.cl <- coeftest(m.vis.cs.wt, vcov. = vcovCL(m.vis.cs.wt, ~ county_state))
m.vis.cs.wt.cl

### Fig 4: Average Marginal Effect Plot----
## Create average marginal effects plot
vis.cme <- plot_cme(model = m.vis.cs.wt,
                    type = "response",
                    variables = "treat",
                    vcov = ~county_state,
                    condition = c("NewGasDistance_s", "post_shock"),
                    draw = FALSE)
vis.cme %>%
  filter(post_shock!=.5) %>%
  mutate(condition2 = ifelse(post_shock==1, "Post-Shale Shock", "Pre-Shale Shock"),
         gas_bin = case_when(
           NewGasDistance_s <= -1 ~ "More Visible",
           NewGasDistance_s > -1 & NewGasDistance_s < 1 ~ "Medium Vis",
           T ~ "Less Visible")) %>%
  group_by(gas_bin, condition2) %>%
  summarize(estimate = mean(estimate),
            conf.low = mean(conf.low),
            conf.high = mean(conf.high))

p.vis <- vis.cme %>%
  filter(post_shock!=.5) %>%
  mutate(condition2 = ifelse(post_shock==1, "Post-Shale Shock", "Pre-Shale Shock")) %>%
  ggplot(aes(x=NewGasDistance_s,y=estimate,fill=condition2,group=condition2)) +
  geom_ribbon(aes(ymin=conf.low,ymax=conf.high),alpha=.8) +
  geom_line() +
  geom_hline(yintercept = 0, color = "grey60") +
  scale_fill_grey() +
  theme_bw(base_size = 14) +
  scale_x_continuous(expand = c(0,0), breaks = seq(-2, 4, 2), limits = c(-2.55,4)) +
  annotate("text", x = .5, y = 10.5, label="Post-Shale Shock", angle = 9.5, size = 7, color = "grey10") +
  annotate("text", x = .5, y = -10, label="Pre-Shale Shock", angle = -4, size = 7, color = "grey60") +
  labs(x="Standardized Distance to New Gas Power Plant",
       y = "Effect on Republican Vote Share (pp)",
       title = "Average Marginal Effect",
       fill = "") +
  theme(
    panel.grid = element_blank(),
    legend.position = "",
    plot.title = element_text(hjust=.5),
    legend.background = element_rect(fill='transparent'),
    legend.box.background = element_rect(fill='transparent',color="transparent"),
    plot.margin = margin(b=-25)
  )
p.vis
p.dist <- subset(g, year >= 1992 & !is.na(treat)) %>%
  ggplot() +
  geom_histogram(aes(x=NewGasDistance_s)) +
  scale_y_continuous(expand = c(0,0), limits = c(0, 750)) +
  scale_x_continuous(expand = c(0,0), breaks = seq(-2, 4, 2), limits = c(-2.55,4)) +
  theme_bw(base_size = 14) +
  theme(panel.grid = element_blank(),
        plot.title = element_text(hjust=.5)
  ) +
  labs(y="", x = "Moderator Distribution", title = "") +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    plot.margin = margin(t=-40,b=-10)
  )
#Combine plot
combined_plot <- p.vis + p.dist +
  plot_spacer() +
  plot_layout(heights = c(3, .6, 0))  # Adjust numbers for desired sizing
combined_plot

## Save Figure 4 output
ggsave(
  device = "tiff",
  filename = here("output", "figures", "fig4_visibilty.tiff"),
  combined_plot,
  dpi = 300,
  width = 5.75,
  height = 4.25,
  scale = 1.5
)
ggsave(
  filename = here("output", "figures", "fig4_visibilty.png"),
  combined_plot,
  dpi = 300,
  width = 5.75,
  height = 4.25,
  scale = 1.5
)

### Table K1----
vis.outcome <-
  g %>%
  filter(!is.na(treat) & year>=1992) %>%
  group_by(fips) %>%
  dplyr::summarize(
    y_avg = mean(RepVotesMajorPercent),
    y_sd = sd(RepVotesMajorPercent)
  ) %>%
  ungroup() %>%
  dplyr::summarize(
    y_avg = round(mean(y_avg),2),
    y_sd = round(mean(y_sd),2)
  )
vis.outcome.s <-
  g %>%
  filter(!is.na(treat) & year>=1992) %>%
  group_by(state) %>%
  dplyr::summarize(
    y_avg = mean(RepVotesMajorPercent),
    y_sd = sd(RepVotesMajorPercent)
  ) %>%
  ungroup() %>%
  dplyr::summarize(
    y_avg = round(mean(y_avg),2),
    y_sd = round(mean(y_sd),2)
  )

vis.coefs <- c(
  "treat"="Coal",
  "post_shock"="Post-Shale Shock",
  "NewGasDistance_s"="Gas Plant Distance",
  "NewGasDistance_c"="Gas Plant Distance",
  "treat:post_shock"="Coal $\\times$ Post-Shale Shock",
  "treat:NewGasDistance_s"="Coal $\\times$ Gas Plant Distance",
  "treat:NewGasDistance_c"="Coal $\\times$ Gas Plant Distance",
  "treat:post_shock:NewGasDistance_s"="Coal $\\times$ Post-Shale Shock $\\times$ Gas Plant Distance",
  "treat:post_shock:NewGasDistance_c"="Coal $\\times$ Post-Shale Shock $\\times$ Gas Plant Distance",
  "oilgas_prop_s"="Hydraulic Fracturing Employment",
  "oilgas_prop"="Hydraulic Fracturing Employment",
  "power_prop"="Power Employment",
  "power_prop_s"="Power Employment",
  "power_prop_c"="Power Employment",
  "unionshare"="Coal Union Share",
  "post_shock:unionshare"="Coal Union Share $\\times$ Post-Shale Shock",
  "white_s"="White",
  "post_shock:white_s"="White $\\times$ Post-Shale Shock",
  "hispanic_s"="Hispanic",
  "post_shock:hispanic_s"="Hispanic $\\times$ Post-Shale Shock",
  "foreign_s"="Foreign-born",
  "post_shock:foreign_s"="Foreign-born $\\times$ Post-Shale Shock",
  "college_s"="College",
  "post_shock:college_s"="College $\\times$ Post-Shale Shock",
  "income99pclog_s"="Income per capita (log)",
  "post_shock:income99pclog_s"="Income per capita (log) $\\times$ Post-Shale Shock",
  "poverty_s"="Poverty",
  "post_shock:poverty_s"="Poverty $\\times$ Post-Shale Shock",
  "rural_nonfarm_s"="Rural",
  "post_shock:rural_nonfarm_s"="Rural $\\times$ Post-Shale Shock",
  "poplog_s"="Population (log)",
  "post_shock:poplog_s"="Population (log) $\\times$ Post-Shale Shock",
  "under40_s"="Under 40 years",
  "post_shock:under40_s"="Under 40 years $\\times$ Post-Shale Shock",
  "work_female_s"="Female workforce",
  "post_shock:work_female_s"="Female workforce $\\times$ Post-Shale Shock",
  "pi" = "Propensity score"
)
file_vis <- here("output", "tables", "si_tab_K1_visibility.txt")
modelsummary(
  list(m.vis.cs, m.vis.cs.power, m.vis.dd, m.vis.dd.pi,m.vis.cs.wt),
  coef_map = vis.coefs,
  vcov = ~ county_state,
  add_rows = data.frame(
    c("Outcome Mean", "Outcome SD", "County Fixed Effects", "State Fixed Effects",
      "Election Fixed Effects", "Covariate Balancing Propensity Score"),
    c(vis.outcome.s$y_avg, vis.outcome.s$y_sd, "No", "Yes", "Yes", "No"),
    c(vis.outcome.s$y_avg, vis.outcome.s$y_sd, "No", "Yes", "Yes", "No"),
    c(vis.outcome$y_avg, vis.outcome$y_sd, "Yes", "No", "Yes", "No"),
    c(vis.outcome$y_avg, vis.outcome$y_sd, "Yes", "No", "Yes", "Yes"),
    c(vis.outcome$y_avg, vis.outcome$y_sd, "No", "Yes", "Yes", "Yes")
  ),
  stars = c("*"=.1,"**"=.05,"***"=.01),
  gof_omit = "IC|RMSE|R2$|Log|Std",
  gof_map = gof_map,
  col.names = c(" ", "(1)","(2)","(3)","(4)", "(5)"),
  output = "latex",
  escape = FALSE,
  fmt = 2
) %>%
  group_rows("Treatment and Moderator:", 1, 12) %>%
  group_rows("Time-Varying Covariates:", 13, 20) %>%
  group_rows("Pretreatment Covariates:", 21, 40) %>%
  cat(., file = file_vis)
fix_txt(file_vis)

## Interaction Effect Diagnostics-----
#subset the data to pre and post
df.post <- subset(g, post_shock == 1 & year>=1992 & !is.na(treat))
df.pre <- subset(g, post_shock == 0 & year>=1992 & !is.na(treat))
df.post <- as.data.frame(df.post)
df.pre <- as.data.frame(df.pre)
### Fig. K3: Linear Interaction Diagnostic Plots ----
#check linearity assumption
raw.pre <-
  interflex(
    estimator = "raw",
    data = df.pre,
    Y = "RepVotesMajorPercent",
    D = "treat",
    X = "NewGasDistance_s",
    Z = c("oilgas_prop", "power_prop", "unionshare"),
    FE = c("state", "year"),
    na.rm = TRUE,
    treat.type = "discrete",
    vcov.type = "cluster",
    cl = "county_state",
    show.grid = FALSE,
    theme.bw = TRUE
  )
raw.pre <- raw.pre +
  labs(x="Standardized Distance to New Gas Power Plant", 
       y = "Republican Presidential Vote Share",
       title = "Pre-Shale Shock") +
  theme_bw(base_size = 14) +
  theme(panel.grid = element_blank())
raw.pre
raw.post <-
  interflex(
    estimator = "raw",
    data = df.post,
    Y = "RepVotesMajorPercent",
    D = "treat",
    X = "NewGasDistance_s",
    Z = c("oilgas_prop", "power_prop", "unionshare"),
    FE = c("state", "year"),
    na.rm = TRUE,
    treat.type = "discrete",
    vcov.type = "cluster",
    cl = "county_state",
    show.grid = FALSE,
    theme.bw = TRUE
  )
raw.post <- raw.post +
  labs(x="Standardized Distance to New Gas Power Plant", 
       y = "Republican Presidential Vote Share",
       title = "Post-Shale Shock") +
  theme_bw(base_size = 14) +
  theme(panel.grid = element_blank())
plot.lin <- gridExtra::grid.arrange(raw.pre, raw.post, ncol = 2)
ggsave(
  plot.lin,
  filename = here("output", "figures", "si_fig_K3_visibility_linearity.png"),
  width = 6.5,
  height = 3.2,
  dpi = 300,
  scale = 1.5
)
### Fig. K4: Binning Estimator ----
bin.pre <-
  interflex(
    estimator = "binning",
    data = df.pre,
    Y = "RepVotesMajorPercent",
    D = "treat",
    X = "NewGasDistance_s",
    Z = c("oilgas_prop", "power_prop", "unionshare"),
    weights = "weights_vis",
    FE = c("state", "year"),
    na.rm = TRUE,
    treat.type = "discrete",
    vcov.type = "cluster",
    cl = "county_state",
    theme.bw = TRUE,
    show.grid = FALSE,
    bin.labs = FALSE,
    xlab = "Standardized Distance to New Gas Power Plant",
    ylab = "Effect on Republican Vote Share",
    main = "Pre-Shale Shock"
  )
bin.pre
bin.post <-
  interflex(
    estimator = "binning",
    data = df.post,
    Y = "RepVotesMajorPercent",
    D = "treat",
    X = "NewGasDistance_s",
    Z = c("oilgas_prop", "power_prop", "unionshare"),
    weights = "weights_vis",
    FE = c("state", "year"),
    na.rm = TRUE,
    treat.type = "discrete",
    vcov.type = "cluster",
    cl = "county_state",
    theme.bw = TRUE,
    show.grid = FALSE,
    bin.labs = FALSE,
    xlab = "Standardized Distance to New Gas Power Plant",
    ylab = "Effect on Republican Vote Share",
    main = "Post-Shale Shock"
  )
bin.out <- gridExtra::grid.arrange(bin.pre$figure, bin.post$figure, ncol = 2)
bin.out
ggsave(
  bin.out,
  filename = here("output", "figures", "si_fig_K4_binning.png"),
  width = 6.5,
  height = 2.8,
  dpi = 300,
  scale = 1.5
)
### Fig. K5: Kernel Estimator-----
kern.pre <-
  interflex(
    estimator = "kernel",
    data = df.pre,
    Y = "RepVotesMajorPercent",
    D = "treat",
    X = "NewGasDistance_c",
    Z = c("oilgas_prop", "power_prop", "unionshare"),
    weights = "weights_vis",
    FE = c("state", "year"),
    na.rm = TRUE,
    treat.type = "discrete",
    vcov.type = "cluster",
    cl = "county_state",
    theme.bw = TRUE,
    show.grid = FALSE,
    bin.labs = FALSE,
    nboots = 2000,
    xlab = "Standardized Distance to New Gas Power Plant",
    ylab = "Effect on Republican Vote Share",
    main = "Pre-Shale Shock"
  )
kern.post <-
  interflex(
    estimator = "kernel",
    data = df.post,
    Y = "RepVotesMajorPercent",
    D = "treat",
    X = "NewGasDistance_c",
    Z = c("oilgas_prop", "power_prop", "unionshare"),
    weights = "weights_vis",
    FE = c("state", "year"),
    na.rm = TRUE,
    treat.type = "discrete",
    vcov.type = "cluster",
    cl = "county_state",
    theme.bw = TRUE,
    show.grid = FALSE,
    bin.labs = FALSE,
    nboots = 2000,
    xlab = "Standardized Distance to New Gas Power Plant",
    ylab = "Effect on Republican Vote Share",
    main = "Post-Shale Shock"
  )
kern.pre
kern.post
kern.out <- gridExtra::grid.arrange(kern.pre$figure, kern.post$figure, ncol = 2)
ggsave(
  kern.out,
  filename = here("output", "figures", "si_fig_K5_kernel.png"),
  width = 6.5,
  height = 2.8,
  dpi = 300,
  scale = 1.5
)

## Fig K1: Within State Variation in Visibility ----
p.vis.measure <- g %>%
  ggplot() +
  geom_histogram(aes(x=NewGasDistance_s, y=after_stat(density))) +
  facet_wrap(~state) +
  theme_bw() +
  labs(x="Standardized Change in Distance",y="Density") +
  theme(panel.grid=element_blank()) +
  scale_y_continuous(expand=c(0,0))
p.vis.measure
ggsave(
  filename = here("output", "figures", "si_fig_K1_visibility_raw.png"),
  p.vis.measure,
  dpi = 300,
  width = 6.5,
  height = 6.5,
  scale = 1.5
)

## Sensitivity Analysis of Interaction Effect ----
#regress the treatment on the covariates
m.gasmod <- lm(NewGasDistance_s ~ treat * post_shock + factor(year) + state +
                 oilgas_prop + power_prop + post_shock * unionshare + 
                 white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s + rural_nonfarm_s + poplog_s + under40_s + work_female_s,
               data = subset(g, year >= 1992 & !is.na(treat))
)
m.gasmod.cl <- coeftest(m.gasmod, vcov. = vcovCL(m.gasmod, ~ county_state))
r2dxj.x <- partial_r2(t_statistic = m.gasmod.cl["poverty_s",3], dof = m.gasmod$df.residual)
r2dxj.x
#regression of the outcome on the treatment and the covariates
r2yxj.dx <- partial_r2(t_statistic = m.vis.cs.power.cl["poverty_s", 3], dof = m.vis.cs.power$df.residual)
r2yxj.dx
#conduct sensitivity analysis
vis.sens <- sensemakr(estimate = m.vis.cs.power.cl["treat:post_shock:NewGasDistance_s",][[1]],
                      se = m.vis.cs.power.cl["treat:post_shock:NewGasDistance_s",][[2]],
                      dof = m.vis.cs.power$df.residual,
                      r2dxj.x = r2dxj.x,
                      r2yxj.dx = r2yxj.dx,
                      kd = c(.43, 2, 4),
                      alpha = .1)
summary(vis.sens)
vis.sens$bounds$bound_label <- paste0(c(".43", "2", "4"), "x Poverty")
### Fig. K6 ----
png(here("output", "figures", "si_fig_K6_sens_vis.png"), width = 6.5, height = 3, units = "in", res = 300, pointsize = 7)
par(mfrow=c(1,2))
plot(vis.sens, sensitivity.of = "t-value", main = "Sensitivity of t-value", lim=.03, lim.y=.03, label.bump.x= 0.003, label.bump.y = 0)
plot(vis.sens, main = "Sensitivity of the Visibility Moderator", lim=.03, lim.y=.03, label.bump.x	= 0.003, label.bump.y = 0)
dev.off()
### Table K2 ----
tabvissens <- here("output", "tables", "si_tab_K2_sense_vis.txt")
vis.sens$sensitivity_stats$treatment <- "Coal $\\times$ Post-Shale Shock $\\times$ Gas Plant Distance"
cat(ovb_minimal_reporting(vis.sens), file=tabvissens)
vis.tab <- readLines(tabvissens)
vis.tab <- vis.tab[-c(1,2,12)]
vis.tab[2] <-  "\\multicolumn{7}{c}{Outcome: Two-Party Republican Presidential Vote Share} \\\\"
vis.tab[3] <- "\\hline"
vis.tab <- c(vis.tab[1], "\\toprule", vis.tab[2:8], "\\bottomrule", vis.tab[9])
cat(vis.tab, file = tabvissens)


## Table K4: Demediation Analysis ----
g.vis <- subset(g, year >= 1992 & !is.na(treat))

## first stage to get the moderating effect of gas distance
first <- lm(
  formula = RepVotesMajorPercent ~ treat * post_shock * NewGasDistance_s + factor(year) + state +
    oilgas_prop + power_prop + post_shock * unionshare +
    white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s + rural_nonfarm_s + poplog_s + under40_s + work_female_s,
  data = g.vis
)

## second stage to get the CDE of visibility net-power plant employment
direct <- lm(
  formula = I(RepVotesMajorPercent - coef(first)["power_prop"]*power_prop - coef(first)["oilgas_prop"]*oilgas_prop) ~ treat * post_shock * NewGasDistance_s + factor(year) + state +
    post_shock * unionshare +
    white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s + rural_nonfarm_s + poplog_s + under40_s + work_female_s,
  data = g.vis
)

# Estimate using sequential-g estimator
g.vis.fips <- subset(g.vis, select = fips) %>% unique()

## block-bootstrap the standard errors
boots <- 1000
fl.boots <- list()
for (b in 1:boots) {
  
  # Get replique
  g.star <- g.vis.fips[sample(1:nrow(g.vis.fips), replace = TRUE), ]
  g.data <- list()
  
  for (j in 1:nrow(g.star)) {
    g.data[[j]] <- subset(g.vis, fips %in% g.star$fips[[j]])
  }
  
  g.data <- do.call(rbind, g.data)
  
  # Estimate first
  boot.first <- lm(
    formula = RepVotesMajorPercent ~ treat * post_shock * NewGasDistance_s + factor(year) + state +
      oilgas_prop + power_prop + post_shock * unionshare +
      white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s + rural_nonfarm_s + poplog_s + under40_s + work_female_s,
    data = g.data
  )
  # Estimate direct
  boot.direct <- lm(
    formula = I(RepVotesMajorPercent - coef(boot.first)["power_prop"]*power_prop - coef(boot.first)["oilgas_prop"]*oilgas_prop) ~ treat * post_shock * NewGasDistance_s + factor(year) + state +
      post_shock * unionshare +
      white_s + hispanic_s + foreign_s + college_s + income99pclog_s + poverty_s + rural_nonfarm_s + poplog_s + under40_s + work_female_s,
    data = g.data
  )
  fl.boots[[b]] <- tidy(boot.direct)
}

# Get data out
boot.out <- bind_rows(fl.boots)

boot.se <- boot.out %>%
  group_by(term) %>%
  summarise(std.error = sd(estimate))
boot.se %>% filter(grepl("treat", term))

first.cl <- coeftest(first, vcov. = vcovCL(first, ~ county_state))

direct.boot <- data.frame(
  term = names(coef(direct)),
  estimate = as.vector(coef(direct))
)
direct.boot <- left_join(direct.boot, boot.se, by = "term")
direct.boot$statistic <- with(direct.boot, estimate / std.error)
direct.boot$p.value <- with(direct.boot, pt(abs(statistic), df = direct$df.residual, lower.tail = FALSE) * 2)
boot.stat <- data.frame("nobs" = nobs(direct), "adj.r.squared" = summary(direct)$adj.r.squared)

boot.tab <- list(tidy = direct.boot, glance = boot.stat)
class(boot.tab) <- "modelsummary_list"

attr(first.cl, "adj.r.squared") <- summary(first)$adj.r.squared

tab_seqg <- here("output", "tables", "si_tab_K4_tab_seqg.txt")
modelsummary(
  list(first, boot.tab),
  stars = c("*"=.1,"**"=.05,"***"=.01),
  gof_map = gof_map,
  gof_omit = "R2|IC|Log|RMSE|Std",
  fmt = 2,
  coef_map = vis.coefs,
  vcov = list(~ county_state, NULL),
  add_rows = data.frame(
    c("State Fixed Effects", "Election Fixed Effects"),
    c("Yes", "Yes"),
    c("Yes", "Yes")
  ),
  escape = FALSE,
  output = "latex"
) %>%
  cat(., file = tab_seqg)
fix_txt(tab_seqg)
